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We consider the isotropic two-dimensional abelian sandpile model from a perspective based 
on two-dimensional (conformal) field theory. We compute lattice correlation functions for various 
cluster variables (at and off criticality) , from which we infer the field-theoretic description in the 
scaling limit. We find a perfect agreement with the predictions of a c = —2 conformal field theory 
and its massive perturbation, thereby providing direct evidence for conformal invariance and more 
generally for a description in terms of a local field theory. The question of the height 2 variable is 
also addressed, with however no definite conclusion yet. 



I. INTRODUCTION 



xysical relevance 



Sandpile models have been invented by Bak, Tang and Wiesenfeld [TJJ as prototypical examples for a class of models 

which show self-organized criticality. The main peculiarity of these models is that they possess a dynamics which 

OO 1 drives them to a critical regime, robust against various perturbations. The ineluctable criticality as well as the 

t-H ' robustness of these specific dynamics could provide a universal explanation of the ubiquity of power laws in natural 

phenomena. Various physical situations have been discussed along this idea, see the recent books j|, 

Sandpile models are among the simplest models showing self-organized criticality. Although their p 

can be questioned, it is believed that they have all the features that should be present in more complicated/physical 

models. Therefore they constitute a useful playground where the most important features can be understood. 

One of the most interesting models is the two-dimensional unoriented abelian sandpile model (ASM) [|l| , which we 

[ — . first briefly recall (recent reviews are P,pl). The model is defined on a L x M square lattice. At each site i, we assign 

a random variable hi, taking its values in the set {1,2,3,4}. We think of hi as a height variable, which counts the 

number of grains of sand at i. Thus a sand configuration is specified by a set of values {hi]i of the height variables. A 

configuration is stable if all hi < 4, and unstable if hi > 4 for one or more sites. The number of stable configurations 

/^ I is equal to 4 LM . 

The discrete dynamics of the model takes a stable configuration C t at time t to another stable configuration Ct+i, 

Q_i! and is defined in two steps. First step is the addition of sand: one grain of sand is dropped on a randomly chosen site 

J ■ of Ct, and this produces a new configuration C' t . The second step is the relaxation to Ct+i- If C' t is stable, we simply 

'-pj set Ct+\ = C' t . If not, the site where hi > 4 topples: it looses 4 grains of sand, and each of its neighbours receives one 

grain, something we write in the form hj — ► hj — Ay for all sites j, with A the discrete laplacian. In the process, one 

neighbour can have its height h > 4, in which case it too topples: it looses 4 grains of sand, each of its neighbours 

receiving one grain. And so on for each site which has a height h > 4, until we reach a stable configuration. Ct+i is 



then set equal to this new stable configuration. The relaxation process is well-defined: it always stops (sand can leave 
the system at the boundaries) and produces the same result Ct+i independently of the order in which the topplings 
are performed (the abelian property). 

One can let an initial distribution over the stable configurations evolve in time according to the dynamics, and 
examine its time limit. Under mild assumptions, one shows H that all initial distributions converge to a well-defined 
and unique distribution P* , call the SOC (for self-organized critical) state. The theory of Markov chains and the 
abelian property allow for a complete characterization of it: P* is uniform on the set 1Z of so-called recurrent 
configurations, and is zero elsewhere (the transient configurations). The number of recurrent configurations is \1Z\ — 
det A ~ (3.21) LM , with A the discrete laplacian on the L x M lattice with open boundary conditions. Although the 
counting of recurrent configurations is easy, the criterion which actually decides whether a given stable configuration 
is recurrent or transient is well-known p,[7| but hard and non-local: in a generic case, one has to scan the whole 
configuration in order to decide whether it is recurrent or not. Explicit calculations are therefore difficult (and few). 

From the point of view of critical systems and conformal field theory, one is interested in the thermodynamic limit 
lini£, ; M^oo P* ■ The result should be a probability measure on the space of spatially unbounded configurations, or 



'Chercheur qualifie FNRS 



equivalently on the infinite collection of random variables hi . Despite the fact that these variables are strongly coupled 
— the couplings are even non-local because of the recurrence condition — , their correlation functions seem to be of 
the usual, local form. In the scaling limit, one could therefore hope to recover a local field theory. 

There are indications that indeed a conformal field theory emerges, like in ordinary critical, equilibrium lattice 
models. In M, a connection with spanning trees was established, which suggests a relationship with the q — limit 
of the g-state Potts model, hence with a c = — 2 conformal field theory, a value confirmed by the calculation of the 
universal finite size correction to the free energy on a finite strip M. The 2-site probability Prob[hi = hj = 1] was 
shown in J8|] to decay algebraically, with an exponent that can be easily accomodated in a c = — 2 free grassmanian 
scalar field theory [ 1. Also the 2-site probabilities for height variables on the boundary of a half-plane domain have 
been computed in [ njfLlJ , and show the same algebraic fall-off as the height 1 variables in the bulk. 

Beyond these concordant elements, no systematic investigation in the sandpile model has been made, to our knowl- 
edge, which can solidly confirm the connection with a c = — 2 conformal field theory. It is our purpose to provide a 
more explicit link between the two. We do this by computing multi-site probabilities of various height variables, and 
by comparing them with the conformal predictions. More specifically, we compute the scaling limit of the 2-, 3- and 
4-site correlations of height 1 variables, but also of other lattice variables, namely finite subconfigurations that can 
be handled by the technique developed in M. 

In fact, we compute these correlations in an off-critical extension of the abelian sandpile model. We evaluate them 
in the scaling regime, extract the scaling limit, and then establish a correspondence with a field theory. In this way 
we strengthen the field-theoretic connection away from criticality, by relating a massive perturbation of the ASM to 
the massive extension of the c = — 2 fcrmionic field theory. One can therefore probe more deeply the structure of both 
pictures, leaving little doubt about the identifications that are to be made. 

The conclusion these calculations allow us to draw is that the c = — 2 theory, and its massive extension, seems 
to provide a field-theoretic description of the height profile of the sandpile model. At least for the cluster variables 
examined in this paper, this is a statement that we could verify explicitly. Other important spatial, non-dynamical 
features of the SOC state must be studied. These include boundary features and avalanche distributions. The latters 
are undoubtedly much more difficult to handle, because they lie at a higher level of non-locality than the height 
variables, since they depend on height values in unbounded regions. Whether they can be accounted for by the 
non-local sectors of the c = — 2 conformal theory remains a largely open question. 

II. LATTICE CALCULATIONS IN THE SANDPILE MODEL 

As recalled above, explicit calculations in the bulk of the lattice are notoriously hard, because of the non-local 
nature of the SOC state (probability measure) P* . 

All four 1-site probabilities Prob[/ii = o], for a = 1,2,3,4, have been computed exactly in the thermodynamic 
limit, but already the calculation for a > 2 |L2| is formidably more complicated than for a — 1 ||. The only 2-site 
probability that has been computed is again for the unit height variables H. 

The technique used to compute the correlation of two unit height variables is a particular case of a beautiful idea 
put forward by Majumdar and Dhar ||. It is based on the important notion of forbidden subconfigurations (FSC), 
and its relation to recurrent configurations. A cluster F of sites, with its heights ft,-, is a FSC if, for each site j E F, 
the number of sites in F and connected to j is bigger or equal to hj. Simple examples of FSCs are two adjacent 
l's (11), a linear arrangement (121), or a cross-shape arrangement with four l's surrounding a central site with any 
height value. A configuration is then recurrent if and only if it contains no FSC |6||7|]. 

The idea used in p) allows to compute the probability of occurrence, in the SOC state, of any cluster that becomes 
an FSC if any of its heights is decreased by one unit. A simple case is a height 2 next to a height 1, but more examples 
are given below in Figure 1. Following [Ol, let us call them weakly allowed cluster variables. 

Let S be such a cluster. The authors oflq] show how one can define a new sandpile model, with its own toppling rules 
(and a new matrix A'), such that the number of its recurrent configurations is the number, in the original model, of 
recurrent configurations which contain S. From this, a simple dcterminantal formula follows, Prob(S') = det A'/ det A. 
Because the new sandpile model is obtained by modifying the original one in the region localized around S, the ratio 
of the two determinants reduces to a finite determinant, even for an infinite lattice. 

This technique has been used to compute the probabilities of various subconfigurations, like those in Figure 1 below. 
The simplest one is the cluster reduced to one site, with height equal to 1. In this case, the new model is obtained 
by changing the toppling rules at 4 sites (the height 1 and three neighbours). A 4-by-4 determinant then yields 
P(l) = Prob[/i, — 1] = -y(l — ^-) ~ 0.074. Allowing for disconnected clusters leads to multi-site correlations such as 



the 2-site correlation of unit heights: 
Prob[hi = hj = 1] = P(l) 2 [l - 



2r 4 



r= |f-i|>l. 



(2.1) 



It was also remarked in || that more general clusters — for instance, a single site with height equal to 2 — can be 
handled using the same ideas, but the corresponding probabilities become infinite series, the terms of which involve 
weakly allowed clusters, of increasing size. Unfortunately, these series seem to be slowly convergent. 
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FIG. 1. On the first two lines are shown the ten smallest weakly allowed cluster variables, up to 
orientations, which contain no more than four sites. Taking the different orientations into account 
makes a total of 57 clusters of weight smaller or equal to 4. In addition, calculations involving the 
four clusters on the last line will be considered in the text. All these clusters will be numbered 
So to S13 from left to right and top to bottom. The reason for including the last two clusters is 



explained in Section VII 



In general, the way the original sandpile model is modified is by removing some of the bonds linking S to its nearest 
neighbourhood, and at the same time, by reducing the threshold at which the sites become unstable (4 in the original 
model), so that the threshold at every site remains equal to its connectivity. These modifications affect all the sites 
of S, plus a certain number of sites which are nearest neighbours of S. All together they form a set we call Ms, the 
cardinal of which depends on the shape of S. The new toppling matrix is then given by A' = A + B, where the 
symmetric matrix B has entries Bij = 1 if the bond linking i to j has been removed, Ba = —n if n bonds off the site 
i have been removed, and is zero otherwise. Then the probability of S (in the original model) is 
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GB) = det (I + GB) 



M s ' 



(2.2) 



Because B is zero outside the finite set Ms, the determinant is finite, in fact of size |Mg|, but requires the knowledge 
of the Green function G = A -1 of the laplacian, at all sites belonging to Ms- 

In the above example where S is just one site with a height equal to 1, the modifications can be pictorially described 
as follows: 

4iL 



4> 



The dashed segments represent the removed bonds, and the numbers on the right lattice indicate the thresholds at 
which the sites become unstable and topple. 

In fact, in the modified lattice shown on the right, the only site to which the 1 is connected has a height bigger 
or equal to 2. So one could as well decrease its height and its threshold by 1, and remove the connection. In this 



way, the site with a height originally equal to 1 is completely cut off from the rest of the lattice, defining a different 
modification of the original ASMQ. Either of them can be used to compute correlations involving heights 1. 

Correspondingly the matrix B that specifies the modifications is a 4-by-4 or a 5-by-5 matrix given by (in an 
obvious ordering) 



B = 




/-3 1 1 1 1 \ 

1-10 ' 
1 0-10 
1 0-10 

V i o o o -1/ 



(2.3) 



For bigger clusters, there is a fair amount of ambiguity in the way the modifications are made in order to freeze the 
cluster heights to what we want. These modifications can affect regions of different sizes, and so can be more or less 
computationally convenient. The least economical solution is the analog of the second modification explained above 
for the unit height. It is also the easiest to describe: one simply cuts the cluster off the rest of the lattice, removing 
all bonds inside the cluster and all connections between the cluster and the outside lattice. There are many other 
choices of intermediate efficiency. For the second cluster in Figure 1 for instance, namely a 2 next to a 1, one may 
consider the following three modifications (among others) : 

I 3 I 3 
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with corresponding B matrices of dimension 8, 7 or 6. For bigger clusters, the difference can be computationally 
noticeable, and so choosing the modifications which affect the smallest possible region makes the calculation of 
determinants easier]]. So one should cut as few links as possible, a prescription which makes sure that the modified 
ASM remains conservative where the original one is: the removal of a bond off one site is always accompanied by the 
lowering by 1 of the threshold at that site, or equivalently the B matrix has row and column sums equal to 0. 

When the cluster S = Ufc S* is disconnected, the matrix B is the direct sum of submatrices B^. The probability 
Prob(S') (the correlation of the subparts Sk) involves the Green function G(i,j) = G(0,i — j) at all sites i,j of S, 
and thus depends on the relative locations and orientations of the various Sk's, and in particular on their separation 
distances. As the original sandpile model is invariant under lattice translation, the probabilities retain the translation 
invariance. For S containing two heights equ al to 1, separated by a distance r, the evaluation of the 8-by-8 determinant 
yields the dominant term r~ 4 given in ( |2.l[ ), independently of the angular distance of the two sites. 

Precisely in the case S is disconnected and contains different pieces separated by large distances, a simple but 
important observation can be made. Because the probability of S is going to depend on the Green function 
G(zk,Zk>) ~ log|zfc — Zk'\ evaluated at points where the subparts are located, one could expect at first sight a 
logarithmic dependence in the separation distances. However due to the property that sand in conserved in the mod- 
ified ASM, the probability in fact only depends on the derivatives (or finite differences) of the Green function. This 
evacuates the logarithmic dependences and turn all correlations into rational series in the various distances {zu — zy)- 

III. THE MASSIVE SANDPILE MODEL 

The previous section summarized the calculation of correlations of cluster variables in the standard ASM. Even 
though it is critical, and self-organized in the dynamical sense, one can drive it off criticality by switching on relevant 



1 Strictly speaking, in this second modified ASM, the removal of the bond connecting the height 1 to its western neighbour 
should be supplemented by the creation of a bond connecting the height 1 to a sink site, so that sand brought by seeding can 
be evacuated. The part of the modifications that affect the sink site plays no role whatsoever, so we may ignore it completely. 
See the second appendix for a detailed argument. 

2 The reader familiar with the technique knows that these determinants can be reduced by appropriate summations of rows 
(or columns). The gain in size is equal to the size of the cluster one considers, but it has its prize, because it renders the entries 
of the reduced determinant more complex. This gain is the same no matter how the ASM is modified. 



perturbations. There are various ways of doing it, but one of the simplest is to add dissipation, whose rate is controlled 
by a parameter t. In effect this introduces a mass m ~ yr, or equivalently a non-zero reduced temperature t. The 
resulting model can be described as a massive (or thermal) perturbation of the massless, critical sandpile model. 
For the purpose of comparing the correlations in the ASM with those of a local field theory, the inclusion of some 
neighbourhood of the critical point is important as it strengthens the connection. 

The way a mass can be introduced in the model is most straightforward, and corresponds to a dissipation of sand 
each time a toppling takes place. We define the perturbed ASM by its toppling matrix (we suppress the explicit 
dependence on x) 

(x Hi = j, 
\,j = < — 1 if i and j are nearest neighbours, (3.1) 

1 otherwise. 

The external driving rate of the sandpile remains the same (one grain per unit of time), but the threshold beyond 
which the sites become unstable is increased from 4 to x. As a consequence, the height variables now take the values 
between 1 and x. Each time a site topples, t — x ~ 4 > grains of sand are dissipated. 

In order to assess the robustness of their SOC features, perturbations of sandpile models have been often discussed, 
with various and sometimes surprising conclusions, see for instance [Oj for a review. Among the many contributions 
on the subject, the reference |15| has been one of the first attempts to see how the non-conservation of sand in the 
top plin g rules can alter the critical properties of the model. In particular, the massive perturbation defined above 



in ( p.l[ ) corresponds to the globally dissipative model studied in that paper, and for which the authors have found 



that the avalanche distributions decay exponentially More recently, the same perturbatio n wa s reconsidered in [TL6l 



in which the exponential decay of the 2-site probability for unit height variables, our Eq. (4.4) below, was proved 



The advantage of the perturbation (3.1) is that it allows the same calculations as the non-perturbed model, in 
the way that has been recalled in Section O. One can in particular compute the correlation functions by the same 
formulae, with however two minor changes. The first one is of course that one uses the massive Green function, with 
a mass fixed by \/x — 4. The second one concerns the B matrices, that define the modified ASM. Because the height 
variables now take the values from 1 up to x, the diagonal entries of B corresponding to sites of the cluster S must 
be set equal to 1 — x, in order to lock the heights into their minimal values. (As a consequence, note that sand is 
conserved at those sites, in the modified model.) 

When doing concrete computations, one needs the value of the Green function at points close to the origin (at sites 
belonging to the same connected subpart), and at points far from the origin (at sites located in different connected 
parts). For the former, one uses a development around t = (in powers of t with logi terms), whereas for the latter, 
one performs a double expansion in inverse powers of the distances, and in (half integral) powers of the perturbing 
parameter t. For arbitrary positions, this development is cumbersome as it depends also on the angular positions. In 
the calculations to be presented in the following sections, we have therefore restricted ourselves to configurations of 
clusters which only require the knowledge of the Green function at points close to a principal or a diagonal axis, for 
which all useful expressions are collected in an appendix. 

The field theory enters as a description of the long distance regime of the ASM correlations (perturbed or not). As 
usual, this requires at the same time an adjustment of the correlation length, or equivalently of the mass. So we are 
interested in computing the scaling regime of correlations. To reach it, we take simultaneously the long distance limit 
R = j l — > oo and the critical limitfj x — 4 = a 2 M 2 — > so that the product ^Jx — AR — > Mr defines the effective mass 
M and the macroscopic distances r. The scale a — > controls the way the limit is taken, and can be thought of as a 
lattice spacing. 

In the actual calculations of correlation functions, large determinants are needed, with entries given by Green 
function values, themselves expressed as power series. In the scaling regime x ~ 4, it is convenient to expand all 
matrix entries and the correlations as power series of \[i. The first non-zero term in a correlation should then be 
directly related to its scaling limit. 

We will finish this section by commenting on the way the calculations have been done, before presenting in the next 
section the results for the unit height variables. 



3 We deliberately take the stand to formally continue all the expressions from integer values of x to arbitrary values x > 4. 
Thus we do not define a family of well-defined sandpile models, parametrized by a real number x > 4. For x rational, this can 
easily be done, however the limit for x going to zero by rational values is not the usual, original model defined for x — 4. We 
susp ect that the model one gets in this specific limit is a model in which the height variables are completely decoupled. See 
[|l7[ for a related discussion. 



Suppose that we want to compute the joined probability for having a certain cluster S at the origin say and an 
other cluster S' at some site i. Each cluster comes with its own set Ms or M$> which contains the sites where the 
ASM has been modified, the modifications themselves being specified by the matrices B and B' . According to the 
discussion of the previous section, this probability is equal to a determinant 

P '° b |S<0).S'(»] = det(,+ (g : go,)(* °))=det(' + G <™ B ,**,). (3-2) 

The G-blocks collectively denotes Green function values evaluated at two sites belonging to the set Ms U Ms' , with 
in addition do = (Goj)*. 

We do not want to know the exact value of this determinant, but rather the terms that are dominant in the scaling 
region, when i is far from the origin. 

Using the standard development of a rank n determinant in terms of the matrix entries, 

det A = ^ e ( (J ) ^1,^(1) ^2,<x(2) ■ • • Ai,a(n) , (3.3) 

o£S„ 



one may distinguish in (3.2) several types of terms. 

The permutations a which do not mix the sites of the cluster 5* with the sites of S' , produce terms which do not 
depend on the distance |i| separating S from S", and thus contribute a term equal to [Prob(S')][Prob(5")]. 

The other permutations necessarily involve an even number of entries from the off-diagonal blocks. As all such 
entries are combinations of Green functions, they decay exponentially with the distance. Therefore the 2-point 
function will be dominated by those terms in the determinant which are quadratic in the off-diagonal Green functions. 
With the help of the formulae in the Appendix A, these Green functions are all reducible to the single G{i) = Go,i> 
and its derivatives. 

The quadratic terms come from the permutations that send one site of the first cluster onto one site of the second 
cluster, and vice-versa (with possibly two other sites). The contributions of all those permutations can be summed 
up to yield a formula written in terms of the minors of the diagonal blocks: 

Prob [5(0), S'(i)] = Prob(S) • Prob(S') - Tr {[Mi(I + GoqB)]* ■ (G 0l B') ■ [Mi(l + GuB'tf ■ (G l0 B)} + ... (3.4) 

Here Mi(j4) = (— l) 4 " 1 ^ det(A~ i -.) denote, up to signs, the minors of A of maximal order (A~ t -■ is the matrix A with the 
i-th row and the j-th column removed). Formula ([Oh is exact modulo quartic, sextic, ... terms in the off-diagonal 
Green functions. It has been used to compute all 2-clustcr correlations considered in this article. 

In order to determine the dominant term in the perturbing parameter i, one still makes an expansion in powers of 
\ft (actually the expansions of elements of the diagonal blocks Goo and Gu involve the two kinds of terms t k > 2 and 
t k / 2 logi). To this end, one develops all Green functions around t = using the formulae of the appendix, and keeps 
the first non-zero term in the trace. Since the mass m or inverse correlation length is related to \/t, a first non-zero 
contribution of the form t( xl+X2 *)/ 2 F (G(iVt), G'(iy/t), . . .)) determines the scaling limit of the correlation, hence the 
corresponding field-theoretic 2-point function, in terms of two fields of scale dimensions x\ and X2, in the usual way. 
In this respect the presence of a logarithmic singularity log t in the final result would be the signal that the scaling 
limit is ill-defined. It turns out, in all the calculations we have performed, that the first non-zero term scales like t 2 
(yielding x\ + x^ = 4). Because the off-diagonal terms start off like \fi — they are differences of Green functions at 
neighbouring sites — , it is enough to expand all Green functions up to order i 3 ' 2 , as has been done in the appendixn 
(the 3-cluster correlations require the expansions to order t 2 ). 

In fact this procedure has anticipated the results on one point. For the purpose of taking the scaling limit, it is 
the dominant term in t that we want to determine, while the above procedure determines the dominant term in t 
among the contributions that are quadratic in the Green functions. So one should also check that no higher than 
quadratic term in the Green functions brings a t 2 contribution. This can easily be done in the following way. Since 
the off-diagonal terms start off like \/t, checking the quartic terms is enough, and one can stop the expansion of the 
off-diagonal blocks to the \fi order. To that order, the two blocks Goi B' and G^o B have all their rows identical. 
Indeed, inside a given column, all entries are finite differences of Green functions evaluated at neighbouring sites, and 



4 It would not be the case if the least economical modifications were chosen (the one that cuts the cluster off the rest of the 
lattice). The B matrix would not have all row sums equal to zero, and consequently the off-diagonal Green functions would 
have non-zero terms of order in t. This would force us to expand everything to order 2 (instead of 3/2) in t. So these 
modifications appear to be doubly inefficient. 



so differ by second order finite differences of Green functions, i.e. by terms of order t. Thus the determinant with 
Goi B' and Gio B as off-diagonal blocks can be reduced to a determinant where the two off-diagonal blocks have but 
their first row non zero, and equal to linear combinations of Green functions. Such a determinant has no term that 
is quartic in the off-diagonal block entries. 



To end this section, we give the expansion analogous to (3.4) that pertains to the calculation of 3-cluster correlations 



Its proof relies on the same arguments as above regarding permutations. For three clusters rooted at sites i,j,k, it 
reads 

Prob[S(i), S'(j), S"{k)] = -2Prob(S) • Prob(S") ■ Prob(S"') 

+ Prob(S) • Prob[S"(j), S"(k)] + Prob(S") • Prob[S(j), S"(k)} + Prob(S"') • Prob[S'(i), S'(j)} 
+ Tr {[Mi(I + GuB)]* ■ (G^B') ■ [Mi(I + G n B')f ■ (G jk B") ■ [Mi(I + G kk B")f ■ (G kl B)} 

+ Tr {[Mi(I + GuB)] 1 ■ {G lk B") ■ [Mi(I + G kk B")f ■ {G kj B') ■ [Mi(I + G^B')] 1 ■ (G Tl B)} + ... (3.5) 

This formula gives all terms of the determinant that are cubic in the off-diagonal Green functions. They are to be 
expanded around t = as discussed above. 

IV. UNIT HEIGHT VARIABLES 

The simplest cluster variable is So, namely the unit height variable. We give in this section its multi-site correlation 
functions, in various configurations, as computed along the lines exposed above. 

The 1-point function, namely the probability a fixed site has height equal to 1, poses no problem (and is in any 
case of little interest for the comparison with a field theory). Making everything very explicit for once, it is given by 



/ /G(0,0) G(1,0) G(l, 0) G(1,0) 

Prob(S'o) = P(l) = dct I - 

V \G(1,0) G(2,0) G(l,l) G(0,0) 
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Here G(m,n) is (A _1 )j for the site i = (m,n), and we have used the symmetries of the Green function. The site 
ordering is O, N, E and S. 

This can be easily computed in terms of complete elliptic functions (see the appendix) , although the result is not 
particularly transparent 

P(l) = ^6 [(* " 4 ) G (°' °) ~ X ] [* 2G ( ' °) ~ 16G ( 1 ' 1 )-( x + 4 )] [i* 2 ~ 8 ) G (°' °) ~ 8G (!' !)-(*- 4 )] 2 • ( 4 - 2 ) 

It goes to 2[2G(1, 1) - 2G(0, 0) + 1] [G(l, 1) - G(0, 0)] 2 = 2(tt - 2)/tt 3 in the limit x -> 4. 

More interesting is its graph, which shows that Prob(So) increases when x goes away from 4 before falling off 
algebraically when x keeps growing. The graph of Prob(So) as a function of x is reproduced in Figure 2 as the 
long-dashed curve. 

A. Two— point correlation 

The joined probability for having a 1 at the origin say and an other 1 at a site i is equal to the 8-by-8 determinant 
Prob[/ lo = / li = l]=det( I + G oo B ,*«*,), (4.3) 



where B is the matrix used in (4.1). Because the two clusters are identical, Gu = Goo- 

As mentioned earlier, the expansion of the Green function at arbitrary points tends to be complicated, so we have 
restricted ourselves to configurations where the Green functions close to a principal axis or a diagonal axis only are 
required. For the 2-site correlation, this only leaves the two possibilities i = (to, 0) and i = (to, m). Using the formula 



(3.4), we found in both cases the same answer: 
Prob[/i = hi = 1] - [P(l)] 2 = 

- t 2 [p(i)] 2 U iC(^|) - \ # (>/*|*|)J^(Vt|»|) + ^ ^o 2 (^KI) + l -^- ffg(V*l*l) \ -- •• • i ^ ' 




FIG. 2. Probabilities of the clusters So up to 64 as functions of the perturbing parameter x. In 
this figure, all probabilities have been normalized to 1 at x — 4 to prevent some of them to melt into 
the horizontal axis. The curves from top to bottom refer to So (long dashes) down to S4 (shortest 
dashes). The two solid lines correspond to S2 and S3. 



with \i\ = m or y/2m depending on whether i is real or on the diagonal. The function Kq is the modified Bessel 
function. Note that the P(l) appearing in the l.h.s. (in the subtraction term) is the off-critical probability, while 
that in the r.h.s. can be taken to be the critical one. 

This formula has a number of instructive and comforting features. The spatial dependence is only through the 
function Kq, that is, the scaling form of the massive lattice Green function. The other functions, denoted Di or Pi in 
Appendix A, and representing the lattice corrections to the scaled, continuum Green function, actually do not enter. 
Moreover the fact that the answer is the same for the two positions of i suggests that the probability is invariant 
under rotations, in agreement with the rotational invariance of the cluster So itself. This is related to the first point, 
since the functions Pi>o and Pi>o represent anisotr opic terms in the lattice Green function. 

Another reassuring feature is that the correlation (4.4) scales like t 2 , to the dominant order, and that all logarithmic 
terms logi have dropped out at that order. This requires massive cancellations because logarithmic terms occur in all 
entries of the bloc ks G qq and Gu , which store Green function values around the origin (see the appendix) . We have 
also checked that ( [4.4| ) is exact up to higher order in t: all terms of order lower than t 2 vanish identically (apart from 
the zero-th order ter m P(l) 2 ), and there are no term quartic or higher in the Green function that would contribute 



a t term. Thus (4.4) is exact to order t 



That the correlation scales like the fou rth power of the mass was expected since the critical correlation decays like 

1-4 



It is easily recovered from (4^4) by taking the limit t —> 



in which the term in K Q 12 



(\i\V~t) 



vw is thc 



only one to survive, reproducing the result (2.1). 

What the above suggests is that the scaled unit height variable goes over, in the scaling limit, to a massive field 0o 
with scale dimension 2, 



lim — ■ 

a^O a 2 



[5{h z/a -l)-P(l)]=<j )0 (z), 
and whose 2-point function reads 



t = a 2 M 2 -> with iy/i = Mz, 



(4.5) 



,, l MV l n [ ^y : = -\r\F { uy\±K "(M\z\)-lK (M\z\)K> Q \M\zV 



2tt 



K' 2 {M\z\) 



1 + 7T 2 

4tt 2 



K 2 {M\z\) .(4.6) 



B. Three— point correlation 



We made the same calculations for the 3-site probability, using the formula ( |3.5| ) • The use of the Green functions 
on the horizontal or the diagonal axis leave essentially two possibilities: either thc three insertion points i, j, and k 
are aligned, or else they form a isoceles right triangle. In both cases, the probabilities scale like t 3 , with all logarithms 
of t cancelled out. The explicit results however differ in these two cases. 



When they form a linear arrangement, be it on the horizontal or diagonal axis, the result for the connected 
probability (i.e. products of lower correlations are subtracted) reads 



Prob[hi = hj 
M 6 



Jaligncd, connected 

[P(l)] 3 x I [K (12) - Kg(12)] [K (13) Kg(13)] [K (23) Kg (23)] + Kg(12) Kg (13) Kg(23) 

Kg(12) K' (13) K' Q (23) - K' Q (12) Kg(13) K' Q (23) + K' (12) K' Q (13) Kg(23) 
K (12) K' (13) K' Q (23) - K' Q {12) K Q (13) K' (23) + K' Q (12) K' (13) K (23) 



- -1 K {12) K (13) K (23) 

■rrO 



(4.7) 



We have written the answer in the scaled form, that is, after the scaling limit in which the sites i,j, k go over to the 
macroscopic positions z\, Z2 and 23. The notation Ka(ij) stands for Ka(M\zi — Zj\). 

For the triangular configuration, we chose the insertion points i = (0, 0) and k — (2m, 0) to be real, and put 
j = (to, to) on the diagonal. The result is slightly different in this case, and reads, in the same notations, 



Prob[hi = hj 
M 6 



Jtriangular, connected 



[P(l)] 3 x ^ 2 Kg (12) K (13) Kg (23) - Kg (12) K (13) K (23) - K (12) K (13) Kg(23) 



Jl 



V2[Kg(12) - K Q (12)} K' (13) K' (23) + V2K' (12) K' (13) [Kg(23) - K (23)} 
+K' (12) [2 Kg (13) -K (13)]K' (23) 

K' (12) K' Q (13) K (23) - K (12) K' (13) K' Q (23)] + -^ K (12) K (13) K Q (23) \. (4.8) 



Exactly the same result was found, as expected, for the rotated configuration where i is at the origin, j = (to, 0) on 
the real axis and k — (m, to) on the diagonal. 

The same comments as for the 2-site correlation apply here but for one point. If indeed the 3-site probability 
scaling ~ £ 3 around the critical point is consistent with the dimension 2 of a unit height variable, one observes that 
the probabilities themselves vanish in the critical limit (M — > 0) . Thus the scaling limit of three unit height variables 
in the usual, unperturbed, ASM vanishes: 



lim Prob [hi = hj = hk = 1] _ 



seal 



2:— 4, connected 



= 0. 



(4.9) 



We have checked this result by using the critical Green functions, and found that the probability for three sites aligned 
along the real axis 



L % 3 fc Jreal, connected 

P(i) 3 r 1 1 1 



~4 ~4 
z 12 z 23 



-12^13 



*23*13 



P(lf 



3P(1) 



f 



1 



1 



-,0 -/O ,-,.£ 

z 12 z 23 z 13 



-,0 -/■£ -/O 

z 12 z 23 z 13 



~£ -,0 -,0 

z 12 z 23 z 13 



4 2 2 
z 12 z 23 z 13 



2 4 2 
z 12 z 23 z 13 



-Yl^Ti^Yi 



higher order (4.10) 



indeed decays like a global power —8 of the separation distances. Moreover the same calculation for the three sites 
aligned on the diagonal axis produces different coefficients. Thus the dominant term of the critical lattice 3-point 
function is not isotropic, contradicting the expected rotational invariance, and so should not survive the scaling limit. 



C. Four— point correlation 

Finally, we have also determined the 4-site probability for unit height variables, at the critical point only, as 
otherwise the number of terms grows quickly. So in this case we have used throughout the calculations the expansions 
at x — 4 of the Green functions, also given in the appendix. 



We have examined two different arrangements of the insertion points, when they are all aligned on the real axis, 
and when they lie at the vertices of a square. 

When they are all aligned on the real axis, the connected 4-site probability takes a very simple form, at the 
dominant order, 



Hi) 4 f 

Prob [hi = kj = h k = hi = 1 real = L±- J 

J connected 4 I 



P(l) 4 f 1 1 1 



(z 12 Z 3 4Z 13 Z 2 4:) 2 (Z13Z 2 4Z 14: Z23) 2 {zi A Z2 3 Z 12 Z 3i ) 



(4.11) 



where the dots represent terms of global power smaller than or equal to —10 (they disappear in the scaling limit), 
and 013 = i — k, . . . (real). 

The other case, for which i = (0,0), j = (m,0), k = (0, m), / = (to, to) are the vertices of a square of side length 
to, is much more rigid as it depends on a single distance m. The result we found for this situation is 

Probfa = hj = h k = hi = 1] «.»» = -- [ —^- + ... (4.12) 

J J connected g m 8 

Before presenting the results for the other cluster variables of Figure 1, we examine the above correlations for the 
unit height random variable from the point of view of the conformal field theory which is the most natural candidate, 
namely the c = — 2 theory, and its massive extension. 

V. CONFORMAL FIELD THEORY 

The c = — 2 conformal field theory has been studied first in the context of polymers Eq], and a bit later served 
as the simplest example of a logarithmic conformal field theory |l9[ . Since then it has been extensively examined by 
many authors ptJ-Eq]. Ref. p4fl in particular presents a clear and rather complete account of the structure the c = — 2 
theory as a rational conformal field theory. Even if it is considered as the simplest situation where logarithms can 
occur, it contains many subtle aspects and probably possesses many different and inequivalent realizations. The one 
which is relevant here is perhaps the most natural one. 

The underlying field theory is formulated in terms of a pair of free grassmanian scalars 9 a = (#, 9) with action 

S=— / e al3 89 a 89 (} = - f 8989, (5.1) 

where e is the canonical symplectic form, e\2 — +1. 

The zero modes of 9, 9, call them £ and £, have been much discussed. Because the action does not depend 
on them, the expectation value of anything that does not contain 9 and 9 explicitely, but only their derivatives 
vanishes identically. In particular the partition function itself vanishes, so the correlation functions are normalized by 
Z' = j T>9'T>9' e~ s , where the primed fields exclude the zero modes £ and £. 

This normalization implies for instance (e a P = — e Q ^) 

(1)=0, (&>=1, (5.2) 

(6 a (z)6 (w))=e a P, (6 a {z)6f 3 (w)U)=e a f } \og\z-w\, (5.3) 

(89 a (z) 89 fj {w)) = , (89 a (z) 89^(w)^) - - — , 2 , (5.4) 

2(z — wy 

(9 a (z 1 )9 (z 2 )9^(z 3 )9\z 4 )) = e a/3 e^ 5 log \z l2 z M \ - e Q V 5 log |zi 3 z 24 | + e aS e^ log \z 14 z 23 \ . (5.5) 

As far as derivatives of fields are concerned — as will be the case in the ASM, at least at the conformal point — , 



one can insert the two zero modes in the correlators, as in (5.4), to take care of the integral on constant fields. The 
functional integral on non-constant fields then yields the usual form for the correlators, obtained from Wick's theorem 
and the kernel of the laplacian. Equivalently, one can define the functional integral for derivative fields by keeping 
the zero modes out, or consider the so-called 7j-£ system Q. 

The stress-energy tensor components T = 2 : 89 89 : and T = 2 : 89 89 : have OPEs characteristic of a conformal 
theory with central charge c = —2. The fields 9 and 9 are primary fields with conformal dimensions (0,0), while the 
bosonic composite field :99: has the following OPE with T 

T(z):99:(w) = — — — ^ + — — -— + •■• (5.6) 



2(z-w) 2 
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It shows that the conformal transformations of : 99 : do not close on itself (and its descendants) but also involve the 
identity and its descendants, which form a conformal module on their own. Thus the identity and : 99 : generate a 
Virasoro module, that is reducible but not fully reducible. This is a characteristic feature of logarithmic conformal 
theories [Q. The field : 99 : is called the logarithmic partner of the identity. It is neither a primary field nor a 
descendant (see below for a field that is primary and descendant without being null). 

The fact that there are two fields with zero scaling dimension is the main source of unusual features (and confusing 
subtleties !), one of them being the existence of two degenerate vacua |0) and |££) (there are two more of fermionic 
nature, |£) and |£)). The above prescription about the insertion of the zero modes can be viewed in the operator 
formalism as the taking of operator matrix elemen ts b etween two distinct in-going and out-going vacua. 



In conclusion, the theory specified by the action (5.1) is a logarithmic conformal theory with central charge c = — 2. 
It contains a non-logarithmic local sector, that retains the central charge value c = —2, and in which derivative fields 
only are considered. Anticipating the analysis to be given below, our results suggest that the ASM scaling fields 
related to height variables precisely lie in this c = — 2 non-logarithmic conformal theory. 

It should also be noted that either theory, logarithmic or non-logarithmic, contains additional non-local (twisted) 
sectors. Although they could play an important role in the sandpilc models, for the description of other lattice 
variables than heights, we will not discuss them here, and refer to |24j for further details. 

We will also need the off-critical, massive extension of the above conformal theory. It corresponds to a perturbation 
by the logarithmic partner of the identity 

1 r M 2 

S(M) = - :d9d9:+ :99:. (5.7) 

7T J 4 

The zero mode problem no longer arises in the massive theory, so that one can normalize the correlation functions 
by the full partition function Z(M) = J ' V9V9 e~ s( - M) . One then obtains 

(9(z)9(w)) = K Q (M\z - w\) , (9(z)9(w)) = (9(z)9(w)) = 0, (5.8) 

and for instance 

M 2 z 
(d9(z) 30(0)) =-—- [Wo(M\z\) - K (M\z\)] . (5.9) 

On account of Kq(x) ~ — log a; for small arguments, the massless limit of the previous equation exists and reproduces 
the expression given in (p^4|) with the zero modes inserted. This is expected since the effect of the zero mode insertion 
is formally to change the normalization factor from Z' to Z. On the other hand, the same does not apply to the 
correlations of the fields 9 a themselves, as the normalizing functional Z(M) goes to zero asAf^O. 

As mentioned above, the cluster variables we consider in this article are all related to derivative fields. The previous 
remark then implies that the off-critical ASM multi-site probabilities have a smooth massless limit, equal to the 
critical probabilities. The scaling form of the off-critical probabilities will be related to the above massive free theory, 
while the critical ones will be computable in terms of the non-logarithmic conformal field theory using the insertion 
prescription. 

VI. SCALING FIELDS FOR CLUSTER VARIABLES 



Let us now reconsider the multi-site probabilities for height 1 computed in the Section IV. The 2-site probability 
suggested that the unit height variable is described by a field with scaling dimension 2, that should be in addition 
scalar since a unit height variable is rotationally invariant. If one assumes that this field is local in 9, 9, the only 
possibilities are :d9 a 89 f3 :, :99 89 a 89' 3 : and M 2 :99:. The second set of fields :99 89 a 89 13 : must be excluded, because 
as explained above, they would produce logarithms in correlation functions, contradicting the observation we made 
in section fil that, in the massless sandpile model (x = 4), the multi-site probabilities are never logarithmic (at least 
those one can compute from the Majumdar-Dhar technique, i.e. from finite determinants). 

It is not difficult to see that 

M 2 
cf>o = -P(l)[: 0989 + 8989: + — :99:\ (6.1) 



is indeed the right combination: its 2-point function is exactly the form given in ( |4.6| ), which was obtained by taking 
the scaling limit of the 2-site probability computed on the lattice. 
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In order to confirm this identification, the field-theoretic 3-point function of 0o can be computed and compared 
with the lattice result. In the same notations as in Section IV, one finds for an arbitrary arrangement of the insertion 
points 

(Mzi) 0o(z 2 ) 0o(^)> = -^r x \ \ ( Z -^- + tc.) K (12) [2# '(13) - K (13)} [2# '(23) - #o(23)] + perm. 

ID Z \Z13Z23 / 



7T-2 



213223 
1-2131 1*23 1 



c.c. ) if (12) if ( 13 ) ^o(23) + perm. 



(6.2) 



( , Z |f 13 ?, 3 , ) «(12) - Xo(12)] # (13) tf (23) + perm. 



if (12)^ (13)^o(23)L 



where the permutations that must be added are the two exchanges z\ «-> 23 and z 2 <-> Z3. One easily checks that 
it repro duce s the 3-site probabilities reported in sec tion [TV] for the two arrangements examined there. The massless 
limit of ( |6.2[ ) vanishes, as it clearly follows from (3.1) for M — 0, since the 3-point function will necessarily involve a 
Wick contraction of a d6 a with some dO 13 . 

Finally the 4-point function can be compared. For convenience, we give the field-theoretic result in the massless 
regime: 



^o(-Zl) 4>0{Z2) 00 (23) <f>o(Z4.)) M=0 = 

F(l) 4 , 



p{iy 



P(l) 



P(l) 



4|2 12 |4|z 34 |4 4|z 13 | 4 |z 24 | 4 4|z 14 |4|z 23 |4 



+ 



{z\2 Z34 Z 13 2 24 ) 2 (Zl3 2 24 Zu Z23) 2 {Zl4 Z 2 3 212 23<l) 2 



C.C. 



(6.3) 



where only the last term within the curly brackets represents the connect ed pa rt of the 4-point function. When the 
four insertions lie on the real axis, it clearly reproduces the lattice result (4.11), and when they are the vertices of a 
square of side m, z\ = 0, z 2 = fn, 23 = im, Z4, = (1 + i)m, it reduces to 



(</>o (21) </>o(2 2 ) ^0(23) 0o(24))m=o, 



, connected 



3 [p(i)Y 



(6.4) 



and again matches the connected 4-site probability ( 4.12J ). 

We believe these comparisons provide enough evidence to assert that the unit height random variable of the sandpile 
model goes over, in the scaling limit, to the field </>o defined in (6.1). In the conformal limit, 0o ~: d0d6+d0d6 := dd :09: 
is a primary field with conformal dimensions (1,1), but is also a descendant of :09:. 

The rest of this section presents analogous results for the other cluster variables pictured in Figure 1. 

We have repeated, for the other thirteen clusters in Figure 1, the same calculations we performed for the unit 
height variable. More precisely, for each of the cluster variables S\ up to S13, we have computed its joint probability 
with a unit height, namely Prob[5o(0), Sk(i)], with i on the principal and on the diagonal axis. From these two 
probabilities one can write down an Ansatz for the field 0^ to which the cluster Sk gets identified in the scaling limit. 
These identifications were subsequently checked to reproduce all 2-site probabilities Prob[Sfc(0), Se(i)], for all pairs 
k,£ = 0, 1, 2, . . . , 13, both on the principal and the diagonal axes. In addition, at least one rotated (or mirrored) 
version of each cluster has been examined, although not systematically (only the correlation with So on both axes). 
The results we found for the rotated clusters are in agreement with the rotations of the fields assigned to the unrotated 
clusters, so that the field of the rotated cluster is the rotated field. Finally, mixed 3-cluster probabilities involving 
unit heights and Si clusters have also been computed. They all confirmed the field identifications. 

All c alculations have been performed exactly, i.e. not numerically. The 2-cluster probabilities take a form similar 
to (4.4), where the coefficients are in general complicated rational expressions of ir. Keeping these coefficients in an 
exact form allows the check of the field identifications to be made in an exact way. For simplicity however, the results 
presented below are given numerically. 

The features of the 2-cluster probabilities are the same as for the unit height variables. We found that all of them 
scale like t 2 , with all logarithmic singularities cancelled out. It implies that all cluster variables go in the scaling limit 
to fields with scaling dimension 2: 



]im±[5(S(i))-P(S)]=M*), i = ~ 

a^o a a 



00 , t = x-4 = a 2 M 2 -> with i\ft = Mi 



(6.5) 



This is somewhat surprising as one could have expected the dimension of the scaling fields to increase with the size 
of the clusters. 
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All cluster variables we have considered have a scaling limit that corresponds to a field of the following form: 

f M 2 - 1 
cj>(z) = -lA:d9d9 + d9d9: + B t :d6d6 + 8986: + iB 2 :d8d9 - 8688: +CP(S) :99:\ . (6.6) 

The (real) coefficients A, Bi, B 2 and C are given in the Table for each cluster. The factor P(S) in front of the term 
: 99 : is the probability of S evaluated at x = 4. Note that the field is not invariant under a rotation of ^ as soon as 
B\ or B 2 is non-zero, but is invariant under a rotation of 7r no matter what the coefficients are. So in particular, the 
scaling limit of the cluster variables does not yield in general conformal fields, but sums of pieces with different tensor 
structures. 



cluster S 


P(S) 


A 


Bi 


B 2 


C 


<$> 


0.0736362 


0.0736362 








1 


<tx±> 


0.0103411 


0.0201433 


-0.00619014 





2 


<2H2><1> 


0.00141994 


0.00449027 


-0.00208908 





3 


<2H2> 


0.00134477 


0.00389417 





-0.000534524 


3 


<2H2X2K^ 


0.00019246 


0.000893234 


-0.000502885 





4 


<2X2H2> 


0.000179829 


0.000752599 


-0.000172832 


-0.000122784 


4 


<|K2> 


0.000173323 


0.000695941 


-0.00012949 


-0.000153467 


4 


-&y-(h- 


0.000179829 


0.000752599 


0.000172832 


-0.000122784 


4 


<2)-(3> 


0.000173106 


0.000692147 


0.000135489 





4 


<|>^<1> 


0.000173106 


0.000692147 


-0.000135489 





4 


<Sio 


0.0000572863 


0.000255127 


-0.0000470539 


-0.0000259802 


5 


Su 


0.00000731457 


0.000042272 


-0.0000130978 


-0.0000071978 


6 


sia= <K><i> 


0.00496687 


0.00969315 


-0.00129393 





2 


Sl3= <±}\> 


0.00404859 


0.010213 





-0.00312009 


2 



TABLE I. For each cluster in Figure 1, the Table gives the values of the parameters A, B 
and C specifying the field that describes the scaling behaviour of the given cluster (see Eq. 
(pTq)). Note in particular that the coefficient C is equal to the size of the cluster. 



As far as numerical values are concerned, the last column of the Table is particularly striking: all entries are integers, 
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simply equal to the size of the cluster. This makes the coefficient of the :89: terms particularly simple and apparently 
regular. The reason for this is unclear. 

The other numbers mentioned in the Table are not in themselves particularly interesting. As mentionned above, 
all these numbers are complicated expressions. For instance the first three numbers on the line corresponding to Sg 
(the last cluster of size 4) are in fact equal to 

2621440 21389312 24279040 14968672 1809776 258037 10061 663 
Prob(S 9 ) - -^ ^- + -^ ^- + -^ ^-^ + — -, (6.7) 

/3tt_8\ / 655360 3389440 2259952 81566 5765 8647 \ 

~ V T 2 J \ 27 7T 5 81 7T 4 + 81 7T 3 9^2" + T^ 96~J' ^ ' ' 

__ /3tt-8\ / 305152 359056 17554 13693 2663 \ 
1 ~\~^~)\8l^ 4 81^" + ~9^ 36^T + _ 96~y ' { } 

A gross feature of the Table is that the (non-zero) numbers are roughly constant for all clusters of the same size, 
namely the probabilities and the coefficients do not change much with the shape of the clusters, but depend essentially 
on their size only. Roughly speaking, these numbers (except C) get divided by 10 when the size increases by one. 

The zeroes in the Table or the equality (up to signs) of coefficients can be understood from the transformations of 
t he clusters and the corresponding fields under the symmetry group of the lattice. One easily sees that the field </> in 
( |6.6[ ) changes under rotations and reflections according to the following rules: 

„s ( (A, — B\,— B%,C) under a ^-rotation, 

(A,B 1 ,B 2 ,C) — ><.)' ' 2 (6.10) 

L (A, B\,—B 2 ,C) under a x- or y-rcflection. 

By convention, all clusters are assumed to be anchored to their lower left site. The rotations are performed about an 
axis passing through that site. 

First of all, the only one to have B\ = B 2 = is the unit height. Indeed it is the only cluster that preserves its 
shape under rotations and reflections, and so one can expect the corresponding field to be a scalar under (continuous) 
rotations and reflections. 

There are clusters whose fields have B 2 = 0, and they are precisely those clusters which are invariant under a 
reflection through the horizontal axis. The same can be said of the rotated clusters for a reflection through the 
vertical axis. That Sg has a coefficient B2 = can be understood along the same lines, although it is not manisfestly 
invariant under reflections. An x-reflection of Sg followed by a rotation by 7r and a translation of two lattice sites bring 
it to itself, except that a height 2 and a height 1 have been swapped. However, the assignment of heights within a 
cluster is irrelevant in the actual computation: the modified ASM is defined in terms of certains bonds being removed. 
Since each site whose height is being constrained looses three out of his four bonds, the actual height assignment is 
irrelevant. In effect, the set Ms which includes all the sites affected by the modifications and the modification matrix 
B itself can be chosen (have been chosen) invariant under a y-reflection. 

In the same way, one sees that S$ and <SV have equal coefficients, up to signs. As represented in the Table and in 
Figure 1, they are related by rotation of ^ and an ^-reflection, with the consequence that their B\ coefficients are 
opposite but the B 2 are equal. The same can be said of Sg and Sg, with the same remark as above regarding the 
locations of the height values within the clusters. 

From these remarks, one easily finds the fields corresponding to different orientations of a cluster. The cluster Sq 
for instance comes in eight different orientations (all anchored to the same site). All of them have the same coefficient 
A ~ 0.000695941 and C = 4, whereas pairs of clusters have coefficients (Bi,B 2 ), or (— B\,B 2 ), or (Bi,— B2) or 
(—B\, —B2). As a consequence, the sum over the corresponding eight fields reduces to a projection onto the scalar 
part, and involves the A and C terms only. 

In a sense, the fact that the fields reflect so well the geometric symmetries of the clusters is surprising. As discussed 
at length in Section O, the actual calculations are based on adequate modifications of the original ASM on a set we 
called Ms, which not only contains sites belonging to the cluster itself, but also sites in its close neighbourhood. Thus 
each cluster drags with itself an invisible shadow, made of the sites in the set Ms \ S. The shadow is a computational 
artefact, but is nevertheless crucial. Moreover it usually breaks or alters the geometric symmetries of the cluster it goes 
with. The insertion of a height 1 for instance, somewhere in the lattice, really requires to consider a 4-cluster pictured 
in Section D. Here the shadow consists of three neighbours of the central site, and clearly breaks the rotational 
invariance. 

We will conclude this section by observing that the height h variables, for h bigger than 4, can be handled in the 
massive ASM exactly like the unit height variables, even more simply. The reason is that a height equal to 5, 6, . . . , x 
can never be in a forbidden subconfiguration, so that the set of recurrent configurations containing a height equal 
to h > 4 at some site i is equal to the set of recurrent configurations on the lattice with i removed. Therefore the 
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modifications needed to freeze the height of a site to h > 4 must simply reduce the threshold at that site to 1, and 
cut it off from the rest of the lattice. This can be implemented by the following matrix 



B 



ll-x 1 1 1 1\ 

' 1 * 

1 

1 

V i o o o 0/ 



(6.11) 



The corresponding probability is simply given by Prob[ft,i = h > 4] = G(0,0), and is logarithmically divergent at 
X = 4. For that reason, one considers instead the probability that hi exceeds 4: 

Prob[h t > 4] = (x - 4) G(0, 0) = fclil K ( -) , (6.12) 

nx \x ) 

which goes to when x — > 4. (The matrix B corresponding to this has —4 as first diagonal entry, rather than 1 — x.) 
if is a complete elliptic function (see Appendix A). 

As for the above clusters, one can compute the correlations of this random variable 8{ht > 4) with itself or with 
the other clusters, and see what field-theoretic description it has in the scaling limit. 

Again the result is simple. The lattice calculation of its own correlation yields 

t 2 
Prob[/i > 4, hi > 4] - Prob[/i > 4] 2 = -— - K$(Vi\i\), (6.13) 

47T 

which suggests the scaling limit 

Z7T 

Correlations with the other cluster variables confirm this limit. It nicely fits the expectation that the field should 
vanish at the critical point. 

VII. THE HEIGHT TWO VARIABLE 

We have so far focused on the class of weakly allowed cluster variables, whose correlations can be handled by 
the technique developped in S, and in turn computed from a finite determinant. The authors point out in that 
article that non weakly allowed cluster variables can in fact be viewed as infinite series of weakly allowed clusters. It 
dramatically complicates their treatment, since a correlation involving a single non weakly allowed cluster requires to 
compute an infinite number of correlations of weakly allowed clusters, of finite but unbounded size. 

In this section, we address the question of the field assignment for the height two variable, in the light of the results 
of the previous sections. We will consider the height two variable, both from the perturbative point of view that we 
have just summarized, and from the conformal point of view. 

That a height two variable can be treated as an infinite sum of weakly allowed cluster variables can be seen as 
follows [pj. Consider the set of recurrent configurations C with a height 2, at the origin say. That set can be divided 
up into two disjoint subsets according to whether the configurations remain recurrent when the 2 is replaced by a 1, 
or become transient upon that replacement. 

The number of those that remain recurrent is the same as the number of recurrent configurations which have a 
height 1 at the origin, because vice-versa, a recurrent configuration with a 1 remains recurrent if the 1 is replaced by 
a 2. So the contribution to P(2) = Prob[/io = 2] from this first subset is exactly equal to P(l). 

For those configurations which become transient, it must be that the 2 belongs to a weakly allowed cluster. This 
weakly allowed cluster can be of various size and shape, and a straight enumeration according to their size leads 
directly to the clusters of Figure 1 (except the first one and the last two) and their various orientations. In this way, 
the second subset is itself divided into a infinite number of disjoint subsets, according to which weakly allowed cluster 
S the height 2 at the origin is part of. The subset labelled by S (fixed size, shape and orientation) contributes P(2) 
a term equal to P(S). 

Putting all together, one obtains, observing that the number P(S) does not depend on the orientation of S, the 
formula 

P(2) = P(1) + J2 P(S) = P(l) + 4,P(S 1 )+AP(S 2 ) + 8P(S 3 )+AP(S i ) + 8P(S 5 ) + ... (7.1) 
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where the summation is over the weakly allowed clusters which are 'anchored' to a height 2. As pointed out in [pj, 
the convergence is very slow. From Table 1, the terms up to Sg furnish the lower bound P(2) > 0.13855, well below 
the exact value P(2) - 0.1739 Jjjj. 

The argument recalled above leading to the perturbative formula for P(2) works similarly for any correlation. The 
result can be expressed as an identity between random variables, 

S(hi-2) = 6(hi-l)+ Y, W))- ( 7 - 2 ) 

w.a.c. S 

Modulo the issue of convergence, this identity is valid when inserted in expectation values. 

The results o f th e previous section suggest that all random variables on the r.h.s. have the same scaling form, given 



by the field in (3.6). Assuming this at all orders, and taking the scaling limit of the previous identity lead to a scaling 



field for the height two of the same form as the scaling field for the height one, namely 

S(hi-2) SC -^5 S a -.8989 + 8989: +PM 2 :99: . (7.3) 

This follows from the observation we made earlier that the other terms 8989 ± 8989 change sign under a rotation by 
~. The sum over the orientations of a cluster make these terms cancel against each other, leaving a scalar field, as it 
should. 

The natural conclusion one could draw from this is that, at the critical point, the heights one and two scale the 
same way and in fact go over, in the scaling limit, to the same — up to normalization — primary field of conformal 
dimensions (1,1). This is the first direct though teneous evidence in favour of such a statement, which has in fact 
been made in |11| , based on an extrapolation to the bulk of a similar statement on the corresponding boundary 
variables, itself relying on the boundary 2-point functions. As plausible and likely as it may be, the extrapolation 
remains uncontrolled, as there are well-known examples of lattice observables that go to different fields, depending 
on whether they lie on a boundary or in the bulk. Thus neither argument is convincing, but both point to the same 
field assignment for the height two variable (and probably similarly for the heights three and four). 

This seems reasonable and likely. It is therefore surprising to observe that it does not appear to be consistent with 
a naive interpretation of the operator product expansions (OPE). To simplify, we consider the critical point, and the 
corresponding conformal field theory. 

The two lattice variables, a height one and a height two, can be taken far apart, and subsequently brought closer 
to each other, until they occupy neighbouring sites, then forming the cluster variable we called Si. In the field- 
theoretic picture, this amounts to taking the two corresponding fields closer and closer to each other, until they 
become coincident, at which point they form a new composite field. The information about what composite fields 
a pair of fields can form when they come close to each other and asymptotically coincident, is contained in their 
operator product expansion (OPE). 

Thus it seems natural to expect that the field assigned to the cluster variable -©-CO" be in the OPE of the field 
corresponding to the height one with the field corresponding to the height two. If one assumes, as argued above, that 
the heights one and two scale to the same field, the required OPE is simply 

: 8989 + 8989 : (z) : 8989 + 8989 : (w) = , * M + '' d f dd '' ^ + :d6d9: W + \ ess sin g U lar, (7.4) 

2\z — w| 4 (z — w) z (z — w) z 

where, fro m di mensional analysis, the less singular terms involve fields of scale dimension strictly larger than 2. One 



sees from (7.4) that the only fields with scale dimension 2 which can be formed in the fusion of a height one with a 
height two are the non-scalar parts of the field making the cluster -@-(l)-. The scalar part of it, : 8989 + 8989 :, is 
missing. (Note that it must be so, since otherwise the unit height variables, represented by that scalar field, would 
have a non-zero (connected) 3-point function.) 

One may observe that the only dimension 2 scalar fields whose fusion produces all field components of -w - ®" are 
logarithmic fields, like -.99(8989 + 8989) :. The change for this logarithmic scalar field has however heavy consequences 
as correlations involving heights two would automatically contain logarithmic functions of the separation distances, 
in addition to the usual rational functions. 

Note that for exactly the same reasons, one could question the field assignment of the height one variable itself, 
despite the fact that the field <j>$ has successfully passed so many tests. Although one cannot bring two heights one 
side by side, one can bring them fairly close to each other, like in the last two clusters of Figure 1 (or Table 1), in 
fact close enough so as not to loose the OPE argument. But then the fields associated to the two clusters S12 and 



S13 must be contained in the fusion of two heights one, i.e. in the fusion (7.4), which we know is not the case. 

Perhaps sandpile models are so special that one would reject the fusion altogether, on the basis that height variables 
have hair ! Because a particular height imposes restrictions on what can stand close to it. For example a height 1 
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forces all its neighbours to be higher or equal to 2, and a height 2 does not allow two of its neighbours to have a 
height 1. This might explain the unconsistency noticed above, but at the same time it denies the very possibility of 
a field assignment. We believe that this issue should be clarified. 

VIII. CONCLUSION 

The power of conformal field theory could bring a much better understanding of the sandpile model, if some of its 
observables could be identified with conformal fields. This is a non-trivial task even for the height variables, which 
are probably the easiest variables to account for in a field-theoretic setting. In addition, and in order to strengthen 
the connection with a field theory, the neighbourhood of the critical point should be investigated. In this article, we 
have taken the first steps towards a systematic study of this relationship, at and off-criticality. 

The off-critical extension of the sandpile that we considered is defined by allowing dissipation, i.e. loss of sand each 
time a site topples. The dissipation rate is controlled by a parameter t > and corresponds to a relevant perturbation 
of the usual abelian undirected sandpile model. 

We have examined multi-site probabilities for the simplest local cluster variables in the off-critical sandpile model. 
By explicit calculations, we have shown that their scaling form can be fully reproduced by a free field theory of massive 
grassmanian scalars. In the massless, critical limit, this theory is a logarithmic conformal field theory with central 
charge c = — 2. The local fields assigned to the various cluster variables however all belong to a non-logarithmic 
bosonic sector. The massive regime, with a mass M ~ \ft directly related to the perturbing parameter, corresponds 
to a thermal perturbation of the conformal theory, i.e. a mass term specified by a logarithmic field. 

We have determined the field assignment for the fourteen cluster variables pictured in Figure 1, and checked their 
consistency against the correlation functions. On the other hand, at the critical point, we have noted a disagreement 
between these assignments and the naive fusion rules of the conformal theory. 

We do not claim that all features of the sandpile models will be comprehensible within a field theory, but some of 
them definitely are. In this respect, other issues than the height variables can be raised: boundary phenomena against 
boundary conformal field theory, the question of the modular non-invariance on a torus (with leaking sites), ... Also 
the relevance and the role of logarithmic fields and twist fields in the c = — 2 logarithmic conformal field theory must 
be further examined. 
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APPENDIX A: GREEN FUNCTIONS 

In this appendix, we collect a number of expressions we have used for the computations of correlations in the 
sandpile model. 

The central object here is the Green function G of the massive discrete laplacian on Z 2 , that is the solution of the 



Poisson equation AG = 1, with A being the finite difference operator given in (3.1). The solution is easily obtained 
by Fourier transform 

G(m > n) = G((m' ) n'),(m + m',n + n')) = — j — , (ro,n),(m / ,n / ) € Z 2 . (Al) 

JJ Att z x — 2 cos k\ — 2 cos k^ 

As explained in the text, values of G are needed at points x which are either close to the origin, or else very far 
from the origin, and in this last case, we have restricted ourselves to points close to a principal or a diagonal axis. We 
treat these three cases in turn. 
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1. The Green function at points close to the origin 



By using the invariance of G under the reflection symmetries of the lattice and its defining equation AG = 1, the 
Green function can be given everywhere in terms of its values on a diagonal. By a suitable change of variables and 
one integration {271 , the diagonal values can be recast into 



(-l) m f* 
G(to, to) = i '— / At 



cos 2mt 



1 - i| sin 2 t 



This can be resolved in terms of the complete elliptic functions 



K(p) 



7T/2 



dt 



1 



Y 1 — p 2 sin t 

E(p) = [ dtJl 
Jo 



q 2 

1 + T 



9g 4 
"64" 



r 7r/2 



— p 2 sin t = 



q 2 13q 4 
~ ~4~ ~ ~64~ 



4 

log{ - 

.9 



21g 4 
128 



3^ 

16 



log( - 
.9 



where q = y/l — p 2 , and where the expansions are given for p < 1 close to 1 . 
In terms of our perturbing parameter £ = x — 4, one finds for instance 

G(0,0) = — #(§), 

7ra; 

o(i,i) = ^{(* 2 -8)^(§)-^(i)}, 

G(2,2) = ^-{(x 4 -16, 2 + 48)^(i)-, 2 (x^8)£;(|)}, 

G(3,3) = -^(( a; 6 -24a ; 4 + 158x 2 -240) X(i) - a; 2 (a; 4 - 16x 2 + 46) £(£)) 
120nx I ' x J 



(A2) 



(A3) 
(A4) 



(A5) 
(A6) 
(A7) 
(A8) 



which can then be expanded around x = 4 by using ( |A3| ) and (A4). They all have the same logarithmic singularity at 
x = 4 as G(0,0), so that the differences G(m,n) — G(0, 0) remain finite when x — > 4. In particular, the critical limit 
of the subtracted diagonal Green function is simply p7J 



lim [G(m,m)-G(0,0)] 

a— >4 



m 1 

7T f-^ 2fc - 1 



(A9) 



fe=i 



2. The Green function on the far diagonal 

For to large, the use of elliptic functions is impractical to extract the asymptotic behaviour in to. Making the 
change of variables z — e lt , the formula (A.2) becomes a integral over a contour that can be deformed to enclose the 



cut lying in between the two roots ±u of the denominator, with u 

,2m 



l(i 



1. This yields 



1 f u 

G(to, to) = / dz 

27r J-u 



v /(z2- U 2 )(z 2_J t) 

The asymptotic expansion of this kind of integral was studied in |29j , from which one finds, using their notations 

1 



(A10) 



G(to, to) 



x(l-u 2 ) y/2~ 7 



. u 2m f 1 | ^> | 3 i2 >~f + , 15 A *> ~ «^ 



5> - e^l> 



3to 64 to 2 



512 



}■ 



where the coefficients A n> are defined from the generating function 



A>{z) 



^/(l + ££*)(! + ££*) 



J2 A ^> zn - 



(All) 



(A12) 
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and are thus themselves infinite (Laurent) series in u 2 , hence in \/i — \Jx — 4. It is not difficult to show that these 
coefficients start off like 



A n > = (-1) 



(2n-l)!l (2 

(2n)!! 



n/2 



o(r 



S-+1 



with the consequence that the to ™ term in (All) takes the form 



(2n-l)l! A n> + 



(-1)" ^( 2n — 1)11 
2" 



(2V2*m) 



[l + scries in i] 



(A13) 



(A14) 



that is, a first term which has the scaling form times corrections in t, independent of the distance m. 
By combining the previous expansion with that of the prefactor of (All), 



1 



1 



e -W2t+mV^V48+... [j + series in j] ^ 



x(l - u 2 ) V27TTO V87rV2iTO 

one eventually finds that the Green function can be written as 

G(m, to) = [l>o(mV2?) + t D 2 (mV2t) + t 2 D 4 (mV2t) + . . . j 



1 mV2P/48+... 



(A15) 



(A16) 



where all functions Di depend on the single scaling vari able m\/2t (the square root of 2 has to be included, since the 
distance from the origin is v2ro). Moreover, from ( A14 ), the first function Do is explicitely given as 



D (z) 



/&KZ 



E 

n=0 



(-1)" ({2n-l)\\ 
2" 



1 



1 



(2z) n 2?r 



ir-Koiz), 



(A17) 



a modified Bessel function. This is to be expected and confirms that the scaling limit of the Green function is indeed 
equal to y~ Ko(Mr), the propagator of a massive scalar. 

For calculations in the ASM model, one still needs the Green functions at points close to the diagonal. The Poisson 
equation is not sufficient, because it would require the knowledge of the Green function all the way down to the 
horizontal axis, but a simple Ansatz similar to (A16|) leads to the following expressions, valid for < k <C to: 



G(to, to + k) 



D Q (z) + kD' (z) x l- + 



D 2 (z) + ^D (z) 



t + 



^D (z) + ^-D' (z)-^D^(z) + ^D' 2 (z] 



2/ :! +...^ e* t/48+ -, (A18) 



where z = m\J2t is the scaled distance. At the order where all the calculations have been performed, the terms shown 
in the previous expression are all that is needed. 



The critical limit of the above expansions is more conveniently computed from (A9) by using the asymptotic 
expansion of the f/'-function pq], or from the integral (A2). The result is 



lim [G(m,m)-G(0,0)] = -— logm- - (£ + log2) - ' 



7 



31 



x^i ' 2n 

with 7 = 0.57721 . . . the Euler constant 



>irm^ 



19207TTO 4 161287TTO 6 



(A19) 



3. The Green function on the far principal axis 



The calculations can be repeated on a principal axis. The integration of (Al) over k 2 followed by the change of 
variable z = e tkl gives G(m, 0) as a contour integral over the unit circle. It can again be deformed to encircle the 
branch cut joining the two roots v < u of the denominator which lie inside the unit circle, yielding 



G(to,0) 



1 

47T 



2,T 



dfei 



^ikirn 



(| — cosfci) — 1 



dz 



(z-u)(z-±)(z-v)(z-l) 



(A20) 
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with u < 1 and v < 1 the two roots of z 2 — (x — 2)z + 1 and z 2 — (x + 2)z + 1 respectively, that is, 



U=-[x-2-y/x(x-4)], V = -[x + 2- y/x(x + 4)]. 



The asymptotic behaviour of this integral for large m can be found again in |2S[| , with the result 

/ ; / ,' f h .-\ i ■ > - - — t , i i ■ > 

G(m,0) = 



47r(l - u 2 ) y/rri 



u'" 1 f 



1 



M-. 



3 A 2> - § 15 A 3> - IA X> 



Am 16 m 2 



04 



(A21) 



(A22) 



The series within the curly brackets is similar to that of the previous subsection, with however -^ substituted for m, 

and with the coefficients A n> as defined in 129], namely by (the coefficients A n> used above correspond to the present 
A n> upon the identification v = —u) 



A>{z) 



1 



J2 A >» zn > 



/(l + i±2» z )(l - i±^z)(l + Mz) 
y v l— tit) ' v 1— v/u ' v 1— u 2 -* 

The usual expansions around ir = 4 now yields 

G(m, 0) = {Po(mVi) + tP 2 (mVt) + t 2 P 4 (mVt) + . . . } e mVF / 24 +- , 

with Po(-z) = 57^0(2:) as before. 

For points close to the horizontal axis, one finds from Poisson equation the expansions for k <C m 



G{m,k) = \P a {z) + 



P2{Z) + -{P,{Z)-P^Z] 



t + 



zt/24+. 



(A23) 



(A24) 



(A25) 



with z = my/t the scaled distance. The ASM calculations also need the values of G(m ± i, k) for small £, and those 
can easily be obtained by expanding the previous result, yielding a Taylor se ries in \ft. 
The critical asymptotic expansion of G(m, 0) can also be computed from (|A20|). One has 



]G(m,0)-G(m + l,0)] 
where F(x) 



l - f * dfci e ik - m F(h) , 
to Jo 



l-e' 



^/(2-cosx) 2 -l ' 

[G(m, 0) - G(m + 1, 0)] 



x =4 Air 

A repeated use of integration by parts then leads to 



1 



x=A 



E 



Ait ^— ' f«m) 
fe>i v y 



(-i) fe r t i 27r 

d^F 



1 f 1 



1 



1 



1 



2tt 1 m 2m 2 2m? 2m 4 



from which one deduces the subtracted Green function itself as 
lim ]G(m,0) -G(0,0)] = -— logm - - (^ + flog 2) 



1 



43 



949 



x— >4 



2tt 
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APPENDIX B: ABOUT THE SINK SITE 



The evacuation of sand is a crucial ingredient to the self-organized criticality of the sandpile models. In order for 
its dynamics to be well-defined — any unstable configuration relaxes to a stable one — , each site should be pathwisc 
connected to a sink, where goes the sand that falls off the pile. The sink is usually omitted in all discussions, perhaps 
because in the ordinary ASM, only the boundary sites are connected to the sink, and that the large volume limit 
takes them to infinity. In the massive ASM however, each site is connected to the sink. One might thus worry about 
a possible role of it in actual computations. 

We show here that the sink has in fact no effect at all and can be omitted completely, be it in the usual or massive 
ASM. The argument is simple and worth being made explicit. 

The discrete dynamics of the ASM recalled in the Introduction uses a toppling matrix that ignores the sink site s. 
To include it, one simply defines an extended toppling matrix A e by adding to A a row and a column, 



1 
V A 



(Bl) 
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The diagonal entry (A ) ss is set equal to 1, in order to freeze the height of the sink site. The rest of the first row is 
equal to 0, since the sink has no connection to the sites of the pile. On the other hand, the first column is not zero: 
Vi >s = —fii if m grains of sand fall off the pile when site i topples. The number rii is equal to n, = — Ylj ^i,j> so 
that all row sums of A c , except the first one, are zero. In the usual ASM, Vi tS is non-zero for the boundary sites only, 
whereas in the massive model, all components are non-zero, with V^ a — 4 — x for all bulk sites. The formula for the 
number of recurrent configurations remains valid, with obviously the same result, dct A c — det A. 

The same method for computing probabilities and correlations of weakly allowed clusters works as before. One uses 
an extended B matrix specifying the way the ASM needs be modified 

B °=(w °b)- ^ 

The first row is clearly always zero, but the first column can be non-zero, depending on the modifications. For those 
called the "least economical" ones in Section O, in which one cuts the cluster off the rest of the lattice, each site i of 
the cluster is left with a sole connection to the sink, so one sets Wi — h — x. 

The probability of a cluster variable S is given by the usual formula, which as before reduces to a finite determinant 

Prob(ff) = det( J Ac ± Bc) = dct(I + A" 1 B c ) = dct (I + A" 1 B c ) . (B3) 

det A c m s u{s} 

However, the restriction to M$ U {s} of A" 1 B c is particularly simple, 

A ^ So= {-A^V A- 1 ) [w b) = {a- 1 W A^bJ' (B4) 

and manifestly leads to the usual result, with no sink, 

Prob(S') = det (I + A" 1 B) . (B5) 
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